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Abstract 

We give a closer look at the Central Limit Theorem (CLT) behavior in quasi- 
stationary states of the Hamiltonian Mean Field model, a paradigmatic one for 
long-range-interacting classical many-body systems. We present new calculations 
which show that, following their time evolution, we can observe and classify three 
kinds of long-standing quasi-stationary states (QSS) with different correlations. The 
frequency of occurrence of each class depends on the size of the system. The dif- 
ferent microsocopic nature of the QSS leads to different dynamical correlations and 
therefore to different results for the observed CLT behavior. 
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1 Introduction 

Very recently there has been a lot of interest in generalizations of the Cen- 
tral Limit Theorem (CLT) [1,2,3] and on their possible (strict or numerically 
approximate) application to systems with long range correlations [4,5], at the 
edge-of-chaos [6], nonlinear dynamical systems the maximal Lyapunov expo- 
nent of which is either exactly zero or tends to vanish in the thermodynamic 
limit (increasingly large systems) [7] , hindering this way mixing and thus the 
application of standard statistical mechanics. A possible application of nonex- 
tensive statistical mechanics [8,9] has been advocated in these cases. Along 
this line we discuss in the present paper a detailed study of a paradigmatic toy 
model for long-range interacting Hamiltonian systems [10,11,12,13,14,15,16], 
i.e. the Hamiltonian Mean Field (HMF) model which has been intensively 
studied in the last years. In a recent article [17], we presented molecular dy- 
namics numerical results for the HMF model showing three kinds of quasi- 
stationary states (QSS) starting from the same water-bag initial condition 
with unitary magnetization (Mq = 1). In the following we present how the 
applicability of the standard or g-generalized CLT is influenced by the dif- 
ferent microscopic dynamics observed in the three classes of QSS. In general, 
averaging over the threee classes can be misleading. Indeed, the frequency of 
appearence of each of these classes depends on the size of the system under 
investigation, and there is no clear evidence that a predominant class exists. 



2 Quasi-stationary behavior in the HMF model 

The HMF model consists of N fully-coupled classical inertial XY spins (rotors) 

'si= {cos 6i, sin 9i) , i = 1, N, with unitary module and mass [10]. One can 
also think of these spins can as rotating particles on the unit circle. The 
Hamiltonian is given by 

N 2 I N 
i — 1 ijj — 1 

where 9i {0 < 9i < 2n) is the angle and Pi the conjugate variable representing 
the rotational velocity of spin i. 

At equilibrium the model can be solved exactly and the solution predicts a 

second order phase transition from a high temperature paramagnetic phase to 
a low temperature ferromagnetic one [10]. The transition occurs at the criti- 
cal temperature T^. = 0.5 which corresponds to a critical energy per particle 
Uc — Ec/N — 0.75. The order parameter is the modulus of the average mag- 
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Fig. 1. Time evolution of the temperature (calculated as twice the average ki- 
netic energy per particle) for three single events representative of the three different 
classes observed at U = 0.69 for initial magnetization Mq = 1. The size of the 
system is iV = 20000. 

netization per spin defined aa M = (1/A^)| E»=i I- Above T„ the spins are 
homogeneously distributed on the circle so that M ~ 0, while below T^, most 
spins are aligned, i.e. rotors are trapped in a single cluster, and M 7^ 0. The 
out-of equilibrium dynamics of the model presents very interesting dynamical 
anomalies. For energy densities U G [0.5,0.75], special classes of initial condi- 
tions such as those called water-bag, characterized by an initial magnetization 
< Mq < 1 and uniform distribution of the momenta, drive the system, after 
a violent relaxation, towards metastable QSS. The latter slowly decay towards 
equilibrium with a lifetime which diverges like a power of the system size N 
[10]. Along the QSS regime, in the microcanonical ensemble, the dynamics 
exhibits a glassy behavior, hierarchical structures, velocity correlations, ag- 
ing, vanishing Lyapunov exponents, etc., and the statistical description of the 
QSS is strongly dependent on the initial conditions [11,12,13]. However, even 
for the same type of initial conditions, we have observed different dynamical 
behaviors [17]. In particular we have found three typical classes of events, as 
illustrated in Fig. 1. Here we plot the temperature of the system (defined as 
twice the average kinetic energy per particle) versus time. We focused on a 
system with size N = 20000, the energy per particle being U = 0.69, for 
which anomalies are more evident. Initial conditions were chosen to be of the 
water-bag kind, with initial magnetization equal to unity (Mq = 1). For details 
about the accuracy of the calculations and the integration algorithm adopted 
see Ref.[12]. Only single events are plotted in Fig. 1, each one representative 
of a given class. One immediately realizes that there are two extreme events 
(which we indicated with 1 and 3) and an intermediate one (indicated as 2) 
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Fig. 2. Relative frequency of occurrence for the three classes of events shown in 
Fig. 1 as a function of A^. A total of m = 20 realizations for each was considered. 
We report also an uncertainty equal to ±A = l/i/m. The three curves add up to 
unity. 

which stays close to event 1 at the beginning, but collapses towards the plateau 
of the event 3 after some time. The relative frequency of occurrence of these 
three classes of events is illustrated in Fig. 2 as a function of the size of the 
system. A total of 20 realizations for each was considered. It is important 
to note that the event of class 3 shown in Fig. 1, and all the events of class 
3 in general, are very similar to the QSS obtained for initial conditions with 
zero initial magnetization (Mq = 0), which are almost homogeneous and have 
very small correlations [12]. In this Lynden-Bell kind of approach, or 

one based on the Vlasov equation has been apphed [15]. Therefore, due to 
the different nature of these QSS, which have in general different microscopic 
correlations, one could expect a different result for the central limit theorem 
behavior shown in Refs. [16,18,19]. 



3 Discussion of numerical results for the CLT 



In this section, following the prescription of the CLT and the procedure 
adopted in Refs. [16,18,19], we construct probability density functions (Pdf) 
of quantities expressed as a finite sum of stochastic variables and we select 
these variables along the deterministics time evolutions of the A^ rotors. More 
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Fig. 3. We present for each class of the QSS found, the different central limit 
theorem behavior observed. A Gaussian (dashed curve) with unitary variance and 
a (7-Gaussian p{x) = Aeq{—(5x'^) with A = 0.66 q = 1.5 and (3 = 1.8 (full curve) 
are also reported for comparison. In the inset a magnification of the central part in 
linear scale is plotted. See text for further details. 

formally, we study the Pdf of the quantity y defined as 
1 " 

y^ = —Y^p^[k) for z = l,2,...,iV , (2) 



where Pi{k), with k = 1,2, ...,n, are the rotational velocities of the ith-iotoi 
taken at fixed intervals of time S along the same trajectory obtained integrating 
the HMF equations of motions. The product S x n gives the total simulation 
time over which the sum of Eq. (2) is calculated. As stressed in Ref. [18], the 
variables y's are also proportional to the time average of the velocities along 
the single rotor trajectories (in fact the 1/y/n scaling is not necessary, and has 
been adopted just to conform to usage). 

In Fig. 3 we show how the shape of the y's Pdf is affected by the existence 
of different classes of events in the QSS of the HMF model. We consider 
the same system of the previous section with = 20000, U = 0.69 and 
Mq = 1 initial conditions. Then we plot the Pdfs corresponding to the three 
different events shown in Fig. 1, the sum (2) being calculated over all the QSS 
extension with S = 100 and n = 2000 (for a total time of 200000 steps, after a 
transient of 200, see orizontal arrow in Fig.l). All the curves were normalized 
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Fig. 4. CLT Pdfs for the same events 1 and 3 of Fig. 3, but with different values of 
6 and n. A Gaussian (dashed curve) with unitary variance and a g-Gaussian with 
A = 0.66, q = 1.5 and P = 1.8 (fuh curve) are also reported for comparison. 

to their standard deviations in order to compare them (notice that, as done 
in [16], in the Eq. (2) we do not subtract an average value calculated over 
the trajectory and all rotors as in [18,19] since it is not important due to this 
kind of normalization). In case of no global correlations, as one would expect 
at equilibrium, the standard CLT holds and the resulting Pdf should be a 
Gaussian (reported as dashed line for comparison). This fact has been verified 
for the equilibrium case in [18]. Actually this is also what happens for the 
event 3 (full circles) and it is most probably due to its similarity with the 
homogeneous QSS recalled in the previous section. On the other hand, if one 
considers the event 1 (gray histogram), strong deviations from the Gaussian 
behavior appears, meaning that global correlations become important. In this 
case a very good fit is obtained with a g-Gaussian defined as 

p{x) = A e,{-f3x^) = [,^(^_i^^^.]V(,-r) > (3) 



where A[q, j3) is a normalizing constant, which coincides in fact with p(0). An 
excellent fitting corresponds to p{y/(j) = 0.66 6^5'^^^^'^^ (full line), a result 
very similar to those discussed in Refs. [18,19]. Finally, the Pdf of event 2 (full 
triangles) lies between the previous two, being (within the time interval that 
has been observed) neither Gaussian nor g-Gaussian but a sort of mixture of 
them, due to the presence of a double QSS plateau in the temperature evo- 
lution of this event. This new scenario, where both Gaussian and g-Gaussian 
attractors appear in the QSS of a given system, implies that, if an average of 
the y^s is performed on many realizations, the particular mixture of the three 
kinds of events discussed here will affect the resulting Pdf, in particular for 
sizes of the system greater than = 10000 (i.e. when the relative frequency 
of events of class 3 becomes important). In this respect, such a scenario also 
qualifies the context within which the results presented in Refs. [18,19] are to 
be considered, which is thus not general, especially for large system sizes. 
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Fig. 5. Different events of class 2 are plotted for the case N=20000. A large vari- 
ability is observed for this class, at variance with the other two. 

As a further remark, in Fig. 4 we test the robustness of the behavior observed 
in Fig. 3 for events 1 and 3. We repeated the calculations on the same events 
of Fig. 3, but with a different time interval 6. In the figure we report the 
cases 5 = 1 (a) and 6 = 4000 (b). One can see that changing 6 and n in 
the calculation of the y's unaffects the Pdf's shape corresponding to those 
events, providing that the total simulation time (i.e. the product 6xn) remains 
unchanged. In this figure, the same Gaussian and g-Gaussian curves of Fig. 3 
are plotted again for comparison. It is important to notice that, for the case 
N = 20000 discussed in this article, we did not study in detail the behavior of 
the CLT for the events of class 2. In this case one could probably distinguish 
a first time evolution, closer to events of class 1, from a final time evolution 
closer to events of class 3. In the present case, however, the first part of the 
plateau is too short to allow the investigation of the CLT behavior. The latter 
can probably be done for systems of larger sizes: we leave this study for the 
future. However, there is another difficulty which is intrinsic with the events 
of this class. In fact they exhibit a general strong variability as shown in Fig. 
5, where we plot the temperature time evolution for several different events 
of class 2 for the case = 20000. Although the qualitative behavior with a 
double plateau is the same, the temperature of the initial plateau can vary. In 
this respect a more detailed investigation is in progress. 

Before concluding, we would like to illustrate that the CLT behavior observed 
for the events of class 3 for Mq = 1 initial conditions is equivalent to that 
observed for any typical event with initial conditions Mq = 0. In Fig. 6 we 
plot the CLT behavior for Mq = in comparison with the case Mq = 1 for 
class 1 events, previously reported in Fig. 3. The Gaussian Pdf obtained for 
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Fig. 6. (a) Comparison of the CLT behavior for the case U = 0.69 initial mag- 
netization Mo = 1 (class 1) vs Mo = 0. The size of the system is = 20000. A 
Gaussian (dashed curve) with unitary variance and a g-Gaussian with A = 0.66, 
q = 1.5 and (3 = 1.8 (full curve) are also reported for comparison, (b) Temperature 
time evolutions of the same events shown in panel (a) . 



the case Mq = is comparable with what obtained in Fig. 3 for the events of 
class 3. 



In Ref. [18], we showed a clear inequivalence between time averages and en- 
semble averages. In Fig. 7, we show that this inequivalence persists also for 
different initial conditions and classes. The ensemble averages of velocities are 
considered at time t = 1000 inside the QSS over 10 events. The Pdfs plotted 
in panels (a) and (b) are to be compared with the time averages curves shown 
in panel (a) of Fig. 6, while those plotted in panels (c) and (d) are to be com- 
pared with the time averages curves of events 1 and 3 of Fig. 3. A Gaussian 
(dashed curve) with unitary variance and a g-Gaussian with A = 0.66, q = 1.5 
and f3 = 1.8 (full curve) are also reported for comparison. It is clear that the 
curves in Fig. 7 are neither Gaussian nor q-Gaussians. 
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Fig. 7. Ensemble averages of velocities for N = 20000, U = 0.69 taken at time 
t = 1000 in the QSS regime for different initial conditions and classes of events. The 
average was done over 10 events. A Gaussian (dashed curve) with unitary variance 
and a g-Gaussian with A = 0.66, q = 1.5 and j3 = 1.8 (full curve) are also reported 
for comparison. See text for further details. 



4 Conclusions 



On the basis of the new calculations presented here, one should distinguish 
among different classes of QSS for a given size of the HMF system. Our tests 
do confirm that correlations can be different for different dynamical realiza- 
tions of the same system, starting from the same class of initial conditions, 
and therefore also the Central Limit Theorem behavior can change. Accord- 
ing to the class considered we can have a Gaussian Pdf, a g-Gaussian one 
or a mixture between the two. In this respect, the presence of a g-Gaussian 
curve confirms the results previously published about the indications for a 
generalized CLT in long-range Hamiltonian systems when ergodicity is broken 
and correlations are strong enough. The present results, where we verify that 
the time evolution of the system is quite sensitive to the specific initial con- 
ditions, reminds what occurs in other anomalous systems, like the Hydrogen 
atom as studied recently in Refs. [20,21]. We have checked also that inequiv- 
alence between time averages and ensemble averages continues to hold within 
the different classes for Mg = 1 and for Mq = initial conditions. Finally, 
although it is true that the frequency of occurrence of events of class 1 and 
2 for very large sizes, and in particular the attractor for the events of class 
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2, remain to be investigated with more accuracy (calculations with better 
statistics are in progress) before advancing definitive conclusions, from these 
simulations one could conjecture that the events of class 1 tend to disappear in 
the thermodynamic limit. However, even if this was the case, their probability 
of occurrence (as well as the probability of occurrence of the first metastable 
plateau of events of class 2) can be significantly different from zero for large 
but finite systems and therefore the applicability of g-statistics in metastable 
states of real complex systems remains a valid and interesting possibility. 
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